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Monte Carlo simulations using the newly proposed Wang-Landau algorithm together with the 
broad histogram relation are performed to study the antiferromagnetic six-state clock model on 
the triangular lattice, which is fully frustrated. We confirm the existence of the magnetic ordering 
belonging to the Kosterlitz-Thouless (KT) type phase transition followed by the chiral ordering which 
occurs at slightly higher temperature. We also observe the lower temperature phase transition of 
KT type due to the discrete symmetry of the clock model. By using finite-size scaling analysis, 
the higher KT temperature T2 and the chiral critical temperature T c are respectively estimated 
as T2 = 0.5154(8) and T c = 0.5194(4). The results are in favor of the double transition scenario. 
The lower KT temperature is estimated as T\ = 0.496(2). Two decay exponents of KT transitions 
corresponding to higher and lower temperatures are respectively estimated as r/2 = 0.25(1) and 
7/1 = 0.13(1), which suggests that the exponents associated with the KT transitions are universal 
even for the frustrated model. 



PACS numbers: 05.70. Jk, 75.10.Hk, 75.40.Mg, 64.60.Fr 



I. INTRODUCTION 

Frustration is one of the interesting subjects in sta- 
tistical physics, mainly because it can induce additional 
symmetry and lead the system to display rich low- 
temperature structures. The so-called two-dimensional 
(2D) fully frustrated XY models have attracted an ex- 
tensive investigation in the last two decades El II II El 
II S II II II 13 IH El 13 III 13 13 113 Due to the 
frustration the systems possess additional discrete reflec- 
tion symmetry Z 2 , apart from the global spin rotation 
symmetry (7(1). The breakdown of these symmetries are 
the onset of two types of phase transitions, namely one 
corresponding to the m agne tic transition of Kosterlitz- 
Thouless (KT) type [3 El and the other to the chiral 
transition. Whether these transitions are decoupled or 
occur at the same temperature has lo ng b een a contro- 
versy I11I10EI1I3I3I3I3I3I3E3 Another 
debated issue is whether the universality class of the chi- 
ral orde ring belongs to the Ising universality class or not 

HHHEai. 

The system has a corresponding physical realization 
on a planar arrays of coupled J osep hson junctions in 
a transverse magnetic field [3 |3j 13 E3 and dis- 
cotic liquid crystals 0. As a 2D frustrated XY sys- 
tem, two lattice systems are frequently studied numeri- 
cally. The first one is the square lattice where the inter- 
actions can be a regular mixture of ferromagnetic (F) 
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and antiferromagnetic (AF) couplings (Villain model) 

nasaslal&ElEiiEi. ^ sec - 

ond one is the AF XY model on the triangular lattice 

H II II El El El El El El E2 

As for the 2D XY model, the effect of the g-fold 
symmetry-breaking fields is an interesting subject |25| : 
that is essentially the same as treating the q-state clock 
model, where only the discrete values are allowed for the 
angle of the XY spins. The U(l) symmetry of the XY 
model is replaced by the discrete C„ symmetry in the 
q-state clock model. It was shown [25] that the 2D q- 
state clock model has two phase transitions of KT type 
at Ti and T 2 (Ti < T 2 ) for q > 4. There is an interme- 
diate XY-like phase between a low-temperature ordered 
phase (T < Ti) and a high-temperature disordered phase 
(T > T2). It is quite interesting to investigate the effect 
of the g-fold symmetry-breaking fields in the case of the 
fully frustrated XY model. Quite recently, Noh et al. 
|H studied the AF six-state clock model on the triangu- 
lar lattice using the Metropolis Monte Carlo simulation 
because of the experimental relevance to CF3Br monolay- 
ers physisorbed on graphite |3|. However, they did not 
pay attention to the lower temperature phase transition 
of KT type. 

It is to be noticed that the existing controversy involves 
very fine values. Most studies claiming single transition 
scenario still do not exclude the possibility of two very 
close critical temperatures. Meanwhile, the studies in fa- 
vor of double transition scenario always found that two 
critical phase transitions occur at slightly different tem- 
peratures. Therefore, it is desirable to obtain precise nu- 
merical information. Recently, much progress has been 
made in the development of efficient algorithms of Monte 
Carlo simulation. Especially, several attempts have been 
proposed for the Monte Carlo algorithms to calculate the 
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energy density of states (DOS ) di rectly. Examples are 
the multicanonical method [28l |29j , the broad histogram 
method j^lj > the flat histogram method 0, H3| , and the 
Wang and Landau method (3j| ■ All of these algorithms 
use the random walk in the energy space. 

In this paper we report our Monte Carlo study on the 
AF six-state clock model on the triangular lattice. The 
ground state (GS) of the AF six-state clock model on the 
triangular lattice has the same structure as the AF XY 
model; therefore this model is regarded as a commensu- 
rate discrete model for the fully frustrated XY model. On 
the other hand, the six-state clock model on the square 
lattice (Villain model) has different GS configurations 
since there exist extra degeneracies. The presence of such 
extra degeneracy may bring about another interest in the 
fully frustrated six-state clock model . However, we will 
not cover such possibility in the present study. The XY 
Villain and the eight-state clock Villain models are com- 
mensurate because they have the same GS configuration. 

For the Monte Carlo method, we employ the Wang- 
Landau algorithm [33j , and the energy DOS is refined by 
the use of the broad histogram relation [34ll35l |. The fact 
that the energy of the six-state clock model is represented 
by the multiple of J/2, where J is the coupling constant, 
is another supporting factor for the study of the six-state 
clock model; it is convenient to treat discrete energy in 
the Monte Carlo simulation of calculating the DOS di- 
rectly. 

The rest of the present paper is organized as follows: In 
the next section we define the model and briefly explain 
the simulation method. Details of the calculation and 
results will be presented in Sec. III. The last section is 
devoted to the concluding remarks. 



II. MODEL AND SIMULATION METHOD 
A. Model and order 

The XY spin model is written with the Hamiltonian 
H = J2 J y ^ ■ S 3 , = J y cos (^ - B o)> (!) 

where (ij) denotes the summation over nearest neigh- 
bor interactions, Si a unit planar spin vector occupying 
the i-th site, and 9i the angle associated with the i-th 
spin. Here, we mainly study the six-state clock model; 
therefore the angle takes discrete values, 9i — 2np/Q with 
p = 0, • • • ,5. The frustration is conveyed by Jij. For the 
Villain model on the square lattice this can be set by 
taking regular mixture of F and AF couplings. For the 
triangular lattice on the other hand, are simply set 
to be uniform AF couplings, Jy = J > 0, so that the 
system becomes fully frustrated. 

The Hamiltonian Q is invariant under the symmetries 
of the global spin rotation U(l) and the global spin reflec- 
tion Zi. The breaking of these symmetries is expected 




FIG. 1: A ground state configuration of the AF six-state 
clock model on the triangular lattice of size 6x6. Spins are 
represented by arrows. Sites belonging to the same sublattice 
have the same orientation of spins. The + and — signs indi- 
cate the handedness of the local chiralities. The ground state 
has 12-fold degeneracy. 

to cause two kinds of ordering, which respectively corre- 
spond to magnetic ordering and chiral ordering. The GS 
configuration is well known as 27r/3-configuration, where 
two neighboring spins align in 2-7r/3 difference in angle, 
which is shown in Fig. ^ We decompose the lattice into 
three interpenetrating sublattices for studying magnetic 
order. A site in a triangle belongs to one of the sublat- 
tices, A, B or C. We assign the magnetic order parameter 
as 

m 2 = (M A 2 + M B 2 + M c 2 ) , (2) 

where Ma = Siga IS the magnetization of sublattice 
A, and N is the number of spins; Mb and Mc follow 
the same definitions for the sublattices B and C. 

To discuss the global spin reflection Z2, we deal with 
the chirality. The local chirality on the elementary tri- 
angle is defined as 

where the z component of the vector chirality is con- 
sidered. The numerical factor in Eq. is chosen such 
that the maximum of the absolute value is one. In the 
GS configuration depicted in Fig. ^ the local chirality 
takes a checkerboard pattern of the right-handed (posi- 
tive) orientation and the left-handed (negative) orienta- 
tion. Then, the staggered chirality 

i 

becomes the order parameter for the Z2 symmetry break- 
ing transition. 

The GS configuration has 12-fold degeneracy which is 
induced by the discrete global spin rotation symmetry 
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Ce with 6-fold and by Z2 symmetry with 2-fold. The 
number of this degeneracy is used as one of the check 
conditions in the calculation of energy DOS. 



B. Simulation method 

We use the Monte Carlo method to calculate the en- 
ergy DOS directly to obtain precise numerical informa- 
tion. First, we briefly describe the Wang-Landau algo- 
rithm 33]. This algorithm is similar to the multicanon- 
ical method (entropic sampling) of Lee [29j . the broad 
histogram method j^lj and the flat histogram method 
[3TI l32| ; the idea is based on the observation that per- 
forming a random walk in energy space with a probability 
proportional to the reciprocal of the DOS, X/g{E), will 
result in a flat histogram of energy distribution. The 
Wang-Landau method introduces a modification factor 
to accelerate the diffusion of the random walk in the early 
stage of the simulation. Since the DOS is not known at 
the beginning, it is simply set g(E) = 1 for all energy E. 
The transition probability from energy E\ to E2 reads 



4000- 



p(Ei — * E2) = min 



g{Ei 
g(E2 



1 



and the DOS g{E) is iteratively updated as 
In g(E) ^ In g(E)+ In fi 



(5) 



(6) 



every time the state is visited. The modification factor 
fi is gradually reduced to unity by checking the 'flatness' 
of the energy histogram; the histogram for all possible E 
is not less than some value of the average histogram, say, 
0.80. 

We also use the broad histogram relation for getting a 
refined DOS. In proposing the broad histogram method, 
Oliveira et al. |30j paid attention to the number of poten- 
tial moves, or the number of the possible energy change, 
N(S, E — ► E'), for a given state S. The total number of 
moves is 



N(S, E —> E + AE) = N 



(7) 



AE 



for a single spin flip process of the Ising model simulation. 
The energy DOS is related to the number of potential 
moves as 

g{E) (N(S, E - E')) E = g(E') (N(S\ E> -> E)) E , , 

(8) 

where (■ ■ ■) E denotes the microcanonical average with 
fixed E. This relation is shown to be valid on general 
grounds 0113) an d we call Eq. © the broad histogram 
relation. We measure the average of the potential move, 
(N(S, E — + E')) E , and use this information for getting 
a better estimate of the energy DOS. It was stressed 
[36l |13 that N(S,E — > E') is a macroscopic quantity, 
which is the advantage of using the number of potential 
moves. We should also note that the broad histogram 
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FIG. 2: Energy DOS of system size L — 48. The energy is 
represented in units of J/2. 



relation does not depend on the particular dynamic rule 
one adopts, and the microcanonical averages of the po- 
tential moves can be obtained by any rule of Monte Carlo 
dynamics. 

In order to reduce calculation time for larger system 
sizes, we break simulation into several energy windows 
and perform random walk in each different range of en- 
ergy. The resultant pieces of the DOS are joined together 
and used to produce the thermal average with the inverse 
temperature p through the standard relation 



(Q) 







J Q(E)g(E)e~P E dE 
jg{E)e-f }E dE ' 



(9) 



Using the parallel machine, we perform the measure- 
ments of the physical quantity Q up to 64 x 10 5 Monte 
Carlo steps. Also, we perform 10 independent runs for 
each system size in order to get better statistics and to 
evaluate statistical errors. 



III. RESULTS 
A. Energy DOS and specific heat 

Here we present the results for the AF six-state clock 
model on the triangular lattice. We have treated the sys- 
tem with the linear sizes L = 24, 36, 48, 60, and 72. We 
apply the periodic boundary conditions. We normalize 
the DOS by using the condition J2e — ® N , and the 
degeneracy in the GS energy, g(Ecs) = 12, is checked 
in order to confirm the accuracy of the calculation. In 
Fig. we show the energy DOS of system size L = 48 
as an example. Here, the energy is represented in units 
of J/2, and the GS energy is given by —(3/2)NJ. 

The energy-dependent data of quantity Q(E) are used 
to calculate the thermal average (Q) by using Eq. @ . We 
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FIG. 3: Temperature dependence of specific heat for system 
size L = 24, 36, 48, 60, and 72. 



calculate the specific heat per spin through the relation 



C(T) 



1 



Nk B T'- 



\{E 2 )-(Ef 



(10) 



where fee is the Boltzmann constant. 

We show the temperature dependence of specific heat 
for different lattice sizes in Fig. |3| The divergent peak 
around T ~ 0.52 in units of J/fce gives a clear sign of 
the existence of second-order phase transition. We also 
observe a hump on the lower temperature side around 
T ~ 0.48, which may be related to the transition of KT 
type. However, we should study the magnetic and chiral 
orders for the detailed analysis of the phase transition. 



B. Correlation ratio 

The critical behavior and the transition temperature 
can be investigated more precisely from the evaluation 
of the order parameter or its corresponding correlation 
function. The magnetic and chirality correlation func- 
tions are defined as the following: 



Gir) = (Si-Si+r), 

7 (r) = (KiKi+r), 



(11) 

(12) 



where r is the fixed distance between spins. Precisely, 
the distance r is a vector, but we have used a simplified 
notation. 

Two of the present authors [38| showed that the ra- 
tio of the correlation functions with different distances 
is a useful estimator for the analysis of the second-order 
phase transition as well as for the KT transition, and this 
correlation ratio can be used for the generalization of the 
probability-changing cluster algorithm |39(. 



At the critical point or on the critical line, the cor- 
relation function g(r) for an infinite system decays as a 
power of r, 



9{r) 



„-(D-2+n) 



(13) 



where D is the spatial dimension and r\ the decay ex- 
ponent. For a finite system in the critical region, the 
correlation function depends on two length ratios, 



g(r,t,L)~r-< D - 2 +^h(r/L,L/0, 



(14) 



where £ is the correlation length. Then, the ratio of the 
correlation functions with different distances has a finite- 
size scaling (FSS) form with a single scaling variable, 



9(r,t,L) 
9(r',t,L) 



f(L/H), 



(15) 



if we fix two ratios, r/L and r/r' . 

In the present work, we set r — L/2 and r' — L/A 
for two distances. Thus, we evaluate the correlation ra- 
tios G{L/2)/G{L/A) and j(L/2)/j(L/A), where G and 
7 are referred respectively to Eqs. l|llf) and i|12|) . It is 
important for two correlated spins to belong to the same 
sublattice; since the fixed distances are set as r = L/2 
and L/A, we choose the system size as a multiple of 12. 



1. Kosterlitz-Thouless transitions 

We show the correlation ratios both for the (a) mag- 
netic and (b) chiral correlations in Fig. 01 From the tem- 
perature dependence of the magnetic correlation ratio 
plotted in Fig. Ufa), we observe that the curves of dif- 
ferent sizes merge in the intermediate temperature range 
(Ti < T < T%), and spray out for the low-temperature 
and high-temperature ranges. This behavior is the same 
as that for the unfrustrated six-state clock model [38| . 
which suggests that there are two phase transitions of 
KT type at T\ and T2. The hump on the lower temper- 
ature side in the specific heat, Fig. |3 may correspond to 
the lower temperature KT transition at T%. The higher 
temperature KT transition at T2 is not obvious from the 
specific heat plot as it is veiled by the divergent peak due 
to the chiral transition. 

We can make a FSS analysis based on the KT form of 
the correlation length, £ cx exp(c/v / t), where t = \T — 
Tkt\/Tkt- The L dependence of Tkt(L) is given by 



Tkt(L) = Tkt + 



(KT 



(ln&L) 2 



(16) 



Using the data of the magnetic correlation ratio R = 
G(L/2)/G(L/A) for different sizes, we estimate two KT 
transition temperatures. We consider the size-dependent 
temperature that gives the constant R. In Fig. EI we plot 
Tkt(L) as a function of l~ 2 with I = InbL for the best- 
fitted parameters in Eq. i|ltj|) . For a fitting function we 
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FIG. 4: Temperature dependence of ratios of the (a) magnetic 
and (b) chiral correlation functions. 

have used a quadratic function in l~ 2 to include correc- 
tion terms. The value of R has been set to be 0.86, 0.88 
and 0.90 for the determination of T2, whereas R has been 
set to be 0.99, 0.985 and 0.98 for T Y . The data with dif- 
ferent R are represented by different marks in Fig.0 but 
they are collapsed on a single curve in this plot, which 
means that b depends on R in Eq. (|16[) and the difference 
of R can be absorbed in the R dependence of b. We es- 
timate the KT temperatures of the magnetic order using 
Eq. O as 

T 2 = 0.5154(8) and 7\ = 0.496(2), 

where the numbers in the parentheses denote the uncer- 
tainty in the last digits. The estimate of T2 is slightly 
lower than the estimate by Noh et al. 26], 0.5175(3). 
It is due to the fact that the moment ratio was used in 
Ref. [26j, and the estimate of the KT temperature be- 
comes higher because of large corrections to FSS |38| . 

Next we consider the decay exponent 77. We first look 
at the constant value of correlation ratio R for different 
sizes and find the associate correlation function G(L/2). 
We give attention to the power-law dependence of the 
correlation function on the system size, G(L/2) ~ L~ v , 
which can be seen from Eq. I|14f> and D is set to be 2. 
We plot G(L/2) versus L for various i?'s in logarithmic 
scale in Fig. |H| The value of 77 is obtained as the slope 
of the best-fitted line for each constant of correlation ra- 
tio. The multiplicative logarithmic corrections for the 
KT transition |19l |40j were shown to be small compared 
to statistical errors. 

We plot 77 thus determined with respect to the fixed 
correlation ratio R in Fig. [7| In the KT phase, R is 
directly related to the temperature. We should note that 
the exponent 77 is meaningful only in the temperature 
range T\ < T < T2 on the fixed line. We show the values 
of R's which give T\ and T2 by arrows in Fig. [7| As can 
be seen, the decay exponent 77 behaves like a typical KT 




FIG. 5: Plot of (a) Ti{L) and (b) Ti(L) of the AF six-state 
clock model on the triangular lattice for L = 24, 36, 48, 60, 
and 72, where I = In bL. The data for H=0.86, 0.88 and 0.90 
are shown by different marks in (a), and those for _R=0.99, 
0.985 and 0.98 in (b). 

transition; that is, the exponent rj continuously changes 
with the temperature in the KT phase. Since rj is almost 
constant for larger R in Fig. [7] the exponent at the lower 
KT temperature T\ is estimated as 

m = 0.13(1). 

For smaller R (higher temperature) side, 77 depends on 
R due to corrections in Fig. Using the fact that the 
fitted value of b in Eq. l|16f) reflects on the difference from 
the transition point, that is, L/£ oc 1/6, we estimate the 
exponent 77 at the higher KT temperature T2 by extrapo- 
lation. The obtained 77 and b are, for example, 0.310 and 
1.76 for R = 0.86, 0.298 and 2.18 for R = 0.88, and 0.284 
and 3.25 for R = 0.90. Plotting the 77's as a function of 
1/6, and extrapolating to X/b — > 0, we obtain 

772 =0.25(1). 

Of course, other dependences such as l/b x are possible; 
such an ambiguity is included in the error. In Fig. [7| 
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FIG. 6: Plot of correlation function G(L/2) versus L. Here 
the slope of the best-fit straight line of each corresponding R 
is the value of exponent 77. 
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tion, 7 (L/2)/ 7 (L/4). 




FIG. 7: Decay exponent r/ of KT phase as a function of mag 
netic coefficient ratio R. Line is just guide to the eyes. 



we show the value of R which gives T2 by the arrow. 
The T) at this R is consistent with the estimated value, 
0.25(1). For the unfrustrated six-state clock model, the 
exponents n> and rji were predicted as 1/4 and 1/9 re- 
spectively ^3 > an d they were confirmed numerically |4l| . 
The present results suggest that the exponents associ- 
ated with the KT transitions are universal even for the 
frustrated model, which the previous work [2(| failed to 
show. 



2. Chiral Phase transition 

The temperature dependence of the chiral correlation 
ratio was also plotted in Fig.^Jb). The existence of chiral 



phase transition can be clearly observed. In the figure, 
there is a single crossing point which indicates the second- 
order phase transition; it corresponds to the divergent 
peak in the specific heat plot. By using the FSS plot 
of chirality correlation ratio, as shown in Fig. |SJ we can 
estimate the critical temperature and exponent v of chiral 
ordering. The estimates are 



T c = 0.5194(4) and v = 0.83(1). 

Our result exhibiting that the chiral transition occurs at 
slightly higher than T2 of KT transition is consistent with 
most studies in favor of double transition scenario. Quite 
recently, Korshunov |42j has discussed that the phase 
transition associated with the unbinding of vortex pairs 
takes place at a lower temperature than the other phase 
transition associated with proliferation of the Ising-type 
domain walls. 

Our estimate for the exp onent v is consistent with the 
results by Lee and Lee |17| and by Ozeki and Ito 0] , but 
contradicts with the result by Olsson ; that is, the crit- 
ical phenomena are not governed by the Ising universality 
class. We have not observed an appreciable size depen- 
dence of the estimated v up to our maximum system size, 
L = 72. Olsson Q argued that corrections to the scaling 
are important in the fully frustrated XY model, and the 
data are consistent with v = 1. Noh et al. also pos- 
tulated that only for large enough system the Ising-like 
behavior is observed. However, using nonequilibrium re- 
laxation study for large enough systems up to L = 2000, 
Ozeki and Ito 01 recently obtained the v = 0.83(2), 
which suggests that the corrections to FSS are not so se- 
rious. Thus, more careful calculations will be needed for 
the critical phenomena of chiral transition. 
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IV. CONCLUDING REMARKS 

In summary, we have investigated the AF six-state 
clock model on the triangular lattice using the Wang- 
Landau method combined with the broad histogram re- 
lation. The model is closely related to the 2D fully frus- 
trated XY model. We have found that the system pos- 
sesses two orderings, spin ordering and chiral ordering. 
The former undergoes the KT transition while the latter 
indicates the second-order transition. We have also ob- 
served the lower temperature KT transition due to the 
discrete symmetry of the clock model. 

Our estimates of the higher KT temperature T2 and 
the critical temperature of chiral ordering T c , that is, 
T 2 = 0.5154(8) and T c = 0.5194(4), support the dou- 
ble transition scenario. The lower KT temperature is 
estimated as T\ = 0.496(2). Two decay exponents 
of KT transitions are estimated as 772 = 0.25(1) and 
771 = 0.13(1), which suggests that the exponents asso- 
ciated with the KT transitions are universal even for the 



frustrated model. 

For the critical phenomena of the chiral transition, our 
estimate of the exponent v, that is, v = 0.83(1), suggests 
that the model does not belong to the Ising universality 
class, but more detailed study is still required. 
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